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where particles interact via soft-shoulder potentials in one dimension. Using a combination of field- 
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liquid by an exotic critical point with central charge c = 3/2. The latter is expression of an 
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the observability of the CLL phase in realistic experimental settings with weakly-dressed Rydberg 
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I. INTRODUCTION 

Frustration plays a fundamental role in our under¬ 
standing of classical statistics mechanics!- 1 One remark¬ 
able example of the effects that frustration can induce on 
a many-body problem is the self-assembly of conglomer¬ 
ate objects, or clusters. The emergence of such compos¬ 
ite objects has been investigated in very different con¬ 
texts, ranging from the physics and chemistry of colloidal 
particles and polymers 21 ^ to two- and three-dimensional 
bosonic systems of ultracold atoms and moleculesP^. 
There, the competition between superfluidity and clus¬ 
tering provides a mechanism to establish supersolidity. 
The latter corresponds to the simultaneous establish¬ 
ment of both diagonal and off-diagonal long-range order, 
a long-sought phenomenon in the context of quantum 
liquid^. 

In the context of one-dimensional (ID) models, liq¬ 
uid phases are usually described by the Tomonaga- 
Luttinger liquid (TLL) universality classP 1 ®. The origin 
of such universality is rooted in the bosonization map¬ 
ping, which allows for the reformulation of interacting 
bosonic, fermionic and spin models onto free bosonic the¬ 
ories. The latter can be exactly solved via path integrals 
as well as conformal field theories technique^-®. TLLs are 
characterized by correlations decaying algebraically as 
a function of distance, implying quasi-long-range-order, 
where the precise values of the exponents depend cru¬ 
cially on the interactions in the microscopic Hamiltonian 
and on its symmetries, such as, e.g., the conservation of 
magnetization or total number of particles. TLLs are 
known to appear in a variety of realizations, including 


edges of topological phased-^, carbon nan otube J^, and 
cold gases of ions 1 ^ atomJis! anc j moleculc d 19 ! 20 ! 

Recently, some of us have proposed that cluster forma¬ 
tion can lead to a new class of quantum liquids, the so- 
called cluster Luttinger liquid (CLL), presenting remark¬ 
ably different features compared to TLlPH One main 
point is that, in the CLL, the essential granularity in the 
liquid is given not by individual particles (as in regular 
TLL), but rather by clusters of particles. As a result, 
in the gapless CLL correlation functions still decay alge¬ 
braically, however, cluster features deriving from the un¬ 
derlying classical cluster structure remain evident. This 
results, for example, in a deformation of the Fermi sur¬ 
face similar to that taking place in Bose metals in ladder 
systemsThis deformation leads to the appearance 
of features in structure factors and the momentum distri¬ 
bution that are not captured by TLLs - where Luttinger 
theorem always holds in the presence of a conserved U(l) 
symmetry associated with particle conservation^®. 

It is the main aim of this work to shed further light onto 
the nature of the CLL state. In particular, we demon¬ 
strate that the CLL and the TLL are distinct phases 
with central charge c = 1 and are separated by an ex¬ 
otic quantum phase transition that displays a quantum 
critical point with c = 3/2. The latter is here associated 
to an emergent supersymmetry, for which we provide nu¬ 
merical evidence. As a model, we focus on hard-core 
bosons trapped in a ID geometry and interacting via a 
soft-shoulder potential. This kind of interaction can be 
treated exactly at the classical level, where it is shown 
to naturally lead to cluster formation. This provides an 
ideal starting point for the analysis of the role of quan- 
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turn fluctuations based on this classical cluster phase. We 
further discuss in detail the possibility to observe these 
peculiar CLL states of matter in experiments with cold 
atomic gases, where soft-shoulder potentials can be engi¬ 
neered by coupling ground- and excited Rydberg states 
with laser light in the weak-dressing regime. There, both 
the range of the soft-shoulder interaction and its height 
can be efficiently controlled via external laser field d -*— 9 - 7 '. 
We conclude by presenting a scheme for the adiabatic 
preparation of the CLL in these systems and an analysis 
of their robustness under typical experimental conditions 
such as finite temperature. 


The paper is organized as follows. In Sec. [TTJ we present 
our model Hamiltonian, which is of the extended Hub¬ 
bard type, and review its exact solution in the classical 
limit. In particular, we first provide a short summary of 
how soft-c ore pot entials can be derived using cold Ryd¬ 
berg atoma^SL 27 and then discuss the param eter r egimes 
where mesoscopic clusters form (see Subsec. IIAI. 


In Sec. |III[ we present a quantum mechanical analysis 
of the strong-cou pling (Subsec. [HI A ) and weak-coupling 
regimes (Subsec. |IIIB| ) of the Hamiltonian. In particu¬ 
lar, using a combination of degenerate perturbation the¬ 
ory from the classical cluster ground state and numer¬ 
ical calculations, in Subsec. Ella] we demonstrate that 
for large enough interaction strengths the effective low- 
energy dynamics is that of a gapless liquid of clusters. In 
Subsec. IIIB we complement this analysis by expanding 
the effective field theory approach presented in Ref. 21 for 
the CLL phase, valid for weak interaction strengths. This 
latter approach retains all cluster constraints exactly and 
thus allows us to make predictions for the correlation 
functions of CLL (see Subsec. Ill CI. Both approaches 
reveal complementary aspects of the physics of cluster 
liquids and their differences with TLLs. 


In Sec. IV we analyze the model numerically using 
the density-matrix-renormalization group (DMRG) al¬ 
gorithm 2829 . The latter is the state-of-the-art method 
to tackle ID systems, and enables the quantitative in¬ 
vestigation of the regime of intermediate interaction 
strengths, where analytical tools are not available. By 
a careful analysis of the scaling of observables such as 
the entanglement entropy (Subse c. |IV A ) and the single¬ 
particle and cluster gaps (Subsec. IVB), we demonstrate 
that the TLL and the CLL are separated by a quantum 
phase transition described by an emergent c = 3/2 con¬ 
formal critical poiifiP^EH] This is an interesting result, 
as the latter points toward the existence of an emergent 
supersymmetry, that is, a supersymmetry that is not ex¬ 
plicitly present in the Hamiltonian. These types of criti¬ 
cal points have a long history in the field of conformal 
field theories and have been recently discussed in the 
context of microscopic spin and fermion ic mo dels with 
many-particle constrainti i 31 l 32 ( In Subsec. IV B 


a numer¬ 


ical analysis of low-lying excitations further demonstrates 
that the spectra of TLLs and CLLs are indeed very dif¬ 
ferent. In fact, while in the TLL the single particle gap 
vanishes in the thermodynamic limit, in the CLL phase 



FIG. 1: Sketch of the degenerate ground state configuration 
for the Hamiltonian in Eq. 0 in the classical regime (t = 
0) at density n = 2/5 [panel (a)] and n = 3/7 [panel (b)]. 
Insets: The structure factor S(q ) versus the lattice quasi¬ 
momentum q for various system sizes N. S(q) is peaked at 
the g-vector characteristic of the formation of ground state 
cluster structures (see main text) and does not shift even for 
very small numbers of particles. 


the gapless degrees of freedom correspond to a vanishing 
cluster gap. Interestingly, we find that the single particle 
gap opens linearly close to the transition point, which we 
discuss to further support the appearance of an Ising de¬ 
gree of freedom in the model, connected to the emergent 
supersymmetry. 

In Secs. |V] and [Vi] we discuss the stability of CLL 
phases in possible experiments with Rydberg-dressed 
atoms-®®!!. In particular, in Sec. [vj we use quan¬ 
tum MonteCarlo (QMC) simulations to demonstrate that 
characteristic features of CLL are evident in correlation 
functions even at relatively high temperatures, signaling 
the stability of the CLL state against thermal fluctua¬ 
tions. In Sec. VI we present a scheme for adiabatic state 
preparation of CLLs by numerically solvi ng t he master 
equation using a quantum jump approach®E3 f or the 
interacting many-body system in the presence of dissipa¬ 
tion. The main conclusion is that dissipation should not 
be detrimental when considering realistic time-scales for 
experiments with weakly-dressed Rydberg-atoms, which 
are usually affected by particle loss due to coupling to 
the Rydberg state. Finally, we draw our conclusions in 
Sec. VII, and present the outlook of our work. 


II. MODEL HAMILTONIAN 

The model Hamiltonian we are interested in describes 
fermionic or hard-core bosonic particles in a ID geometry 
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and reads 

_ r c 

H = i+h.c.) + V EE niUi+t. (1) 

i i i=\ 

Here, b\ ( bi ) are bosonic/fermionic creation (annihilation) 
operators at the site i, nt = b\bi , and t is the tunneling 
rate. 

Interactions between particles have a soft-core profile, 
with depth V and radius r c . These particular features of 
the interaction potential can be experimentally realised 
with Rydberg-dressed cold gases. In the weak-dressing 
regime -C |A|, atoms in their ground states are off- 
resonantly coupled with Rabi frequency fl and red detun¬ 
ing A < 0 to an high-lying Rydberg state. The resulting 
effective potential as a function of the relative distance x 
between pair of atoms read^ 


Q 4 r 6 

V{x) = —y - c ■ 

8A 3 r® + a; 6 


( 2 ) 


where r c = [C§/{2\ A|)] 1 / 6 is the Condon radius and Cq 
is the van-der-Waals coefficient of the addressed Ryd¬ 
berg state. At large distances x r c , V(x) reduces 
to the usual repulsive van-der-Waals interaction between 
Rydberg atoms oc x~ 6 , suppressed by a factor [Q/(2A)] 4 , 
since only a small fraction [ST/(2 A)] 2 of the Rydberg state 
is admixed into the original ground state. However, for 
x < r c , a double Rydberg excitation is prevented by the 
dipole blockade and V(x) saturates to a universal con¬ 
stant value H 4 /8A 3 . 


A. Cluster exchange model 


However, the ground state is not a unique configuration 
but there are exponentially many degenerate configura¬ 
tions with the same energy. These configurations can 
be represented by dividing the system in blocks made 
of particles and holes, where the number of holes is the 
same in each block. Here, e.g., blocks labeled A consist 
of two particles and two holes with a total length of 4, 
whereas blocks labeled B consist of one particle and two 
holes with length 3. In this block model, each density of 
particles n corresponds to a ratio of number of blocks A 
(B), Na(N b ). For example, n = 2/5 corresponds to a 
ratio Na/Nb = 1/2 and n = 3/7 to Na/N b = 1, respec¬ 
tively. The ground state then consists of all permutations 
of blocks A and B (e.g. for n = 3/7 typical configura¬ 
tions are [AABAAB...], [ABAAAB...], ...jAfl. The asso¬ 
ciated ground state degeneracy, assuming open bound¬ 
ary conditions, is d = M\/[(M/ 3)!(2M/3)!] for n = 3/7 
and d = M\/[(M/ 2)!] 2 for n = 2/5, respectively, where 
M = (L — N)/r c is the total number of clusters in the 
system. 


III. THE CLUSTER LUTTINGER LIQUID 
PHASE 

The complex structure of the classical ground states re¬ 
sults in the emergence of an exotic quantum liquid once 
quantum fluctuations are introduced. In this Section, we 
will first present a strong coupling approach in the b >f 
limit, and then a modified bosonization treatment which 
embodies the cluster constraints derived in the classical 
limit. The combination of the two approaches allows us 
to gain a qualitative analytical understanding of the full 
phase diagram of Eq. |l]) , which we quantitatively inves¬ 
tigate in Sec. IV below. 


Before presenting results for the quantum phases of 
Hamiltonian Eq. (jT]), we first summarize results for the 
classical ground state. Particles with soft-shoulder in¬ 
teractions and their cluster-phases have been extensively 
studied in the classical regime in the context of soft mat¬ 
ter physics 2 3 . In one dimensional chains, we study the 
classical ground state from applying the cluster exchange 
model introduced in Ref. m, The two relevant length 
scales in this regime are the cut-off radius r c and the av¬ 
erage distance between particles r* = 1/n, with n = N/L 
the particle density. Here, N and L are the number of 
particles and of lattice sites respectively. This leads to 
three possible regimes in the classical limit: (i) liquid 
r* > r c: (ii) crystal r* = r c and (iii) cluster liquid phase 
r* < r c . In this work we will consider r c = 2 and den¬ 
sities n = 2/5 and n = 3/7 which both correspond to a 
cluster phase. 

The cluster exchange model is illustrated in Fig. |Tj 
For, e.g., r c = 2 and n = 2/5 the state with lowest 
energy consists of single particles and clusters of size 2 
[depicted in Fig. [lja)] with a total energy of E 0 = Vn. 


A. Strong-coupling approach to cluster manifolds: 
the XX model 

Once the degenerate manifold of clustered ground 
states has been identified, it is possible to derive an ef¬ 
fective Hamiltonian describing the dynamics in the limit 
t <C V. In order to do so, we define as V the projector 
operator on the classical ground state manifold, and ap¬ 
ply conventional second order perturbation theory (odd 
order corrections are not present for the case r c = 2) 

H eS . -Ho + VHtQ—*— QH t V + 0(t 4 /V 3 ), (3) 

e - Ho 

where Q = 1 — V, and e is the classical ground state 
energy within the cluster manifold. First, we define ef¬ 
fective spin-1/2 operators Sj as follows: Ordering the 
cluster configuration with an index j £ [1,M], we as¬ 
sociate to the position of each A-type cluster a spin-up, 
and for each R-type, a spin-down. This is a one-to-one 
mapping of the Hilbert space defined by V to the Hilbert 
space of a spin-1/2 chain with M sites. The effective 
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FIG. 2: Main Panel. Energy of the ground state (GS) as a 
function of V/t in the perturbative limit t<F (here t = 1). 
Symbols denote the exact GS energies obtained via exact di- 
agonalization of Eq. 0 . while lines denote perturbative es¬ 
timates for different cluster configurations and system sizes. 
Inset: Absolute value of the energy difference e between the 
exact and the perturbative estimate. The energies are quite 
accurate down to low values of V/t < 10, where the strong 
coupling expansion around the cluster manifold first starts to 
break down. The thin, black line is a guide for the eye in¬ 
dicating a power law decay t A /V 3 . The latter demonstrates 
that higher order terms are irrelevant at large couplings, but 
play a role around V/t ~ 10. The small kink around V/t ~ 2 
is due to the fact that e changes sign in this region. 


Hamiltonian can be then cast in a compact form as a 
spin chain. For the r c = 2 case, diagonal contributions 
of the type Sf Sf +1 do not contribute at lowest order and 
we get 

,2 m 

H eS . Y,K S t S I+i + h - c ~) + 2 ]- ( 4 ) 

i=i 

The strong coupling limit can then be mapped to a sys¬ 
tem of hard-core bosons (spin-1/2) hopping in an arti¬ 
ficial lattice created by the underlying cluster structure. 
This confirms that the CLL, which is adiabatically con¬ 
nected to the strong coupling limit, is indeed described 
at low-energies by a c = 1 conformal field theory (CFT), 
which is interpreted as a Luttinger liquid of composite 
cluster particle^. 

In order to benchmark the validity of the strong cou¬ 
pling expansion and to check to which extent correc¬ 
tions to Eq. Q are quantitatively relevant, we compare 
in Fig. [2] the energy of a system of L = 14 and 10,20 
sites for n = 3/7 and 2/5, respectively, and r c = 2, us¬ 
ing exact diagonalizations of the full Hamiltonian as well 
as perturbative estimates. We find a very good agree¬ 
ment between the exact and perturbative results down 
to relatively small values of V/t ~ 5. This indicates that 
all relevant quantum dynamics is well described within 


our strong-coupling model detailed above, where clusters 
play the role of the mesoscopic degrees of freedom. At 
smaller interaction values, higher order terms are non- 
negligible, and moreover the cluster assumption breaks 
down, as demonstrated by the poor agreement between 
the exact and perturbative results below V/t ~ 5 [see 
inset in Fig. [2]. 


B. Beyond perturbation theory: Low-energy field 
theory of cluster Luttinger liquids 

While perturbation theory provides an understanding 
of the large V limit, an analytical picture at intermediate 
couplings is hindered by the complex structure of the 
interactions. In Ref. EH we discussed how cluster-type 
constraints can be generically applied to Haldane’s 
bosonization prescription to derive a modified mapping 
between the original microscopies fields and continuous 
bosonic variables. The main feature of this treatment 
is that it captures the correct behavior of correlation 
functions and the deformation of the Bose surface, 
which is instead not accessible by direct bosonization of 
Eq. ([I]). Here, we provide a detailed derivation of the 
mapping and a discussion of its consequences on various 
observables. Analytical predictions will be compared 
with exact numerical simulations in the next sections. 

We start from the second step of Haldane’s construc- 
tiorP2l, by considering a more complicated shape of the 
initial particle distribution that can reflect the underly¬ 
ing cluster structure. After taking the continuum limit, 
this can be cast as the following constraint on the density 
distribution for the density operator p 

N M 

p(x ) = ^2 5(x - X n ) = ^2 5(x - X m ) + ^2 K x ~ x l) = 

n =1 m =1 ££cl 

M 

= ^2 f(x m )S( x - Xm)- (5) 

m= 1 

Here, Yl x f( x m ) = N, and the function f(x m ) G {1, 2} 
describes both one- and two-particle cluster structures. 

In the center equality, ^2 cl represents the sum over all 
two-particle clusters. Formally, the summation Y2 c i con “ 
straints the values of xe such that there exists an m such 
that Xu = x m W. The key point here is that the sum is 
split into two parts: the first one describes the fact that 
there are some clusters where at least one particle resides, 
while the second sum takes into account the possibility 
of having clusters formed by two particles. We now in¬ 
troduce a new field, tp c i(x), which accounts for quantum 
fluctuations of the cluster density. The field satisfies: 

<Pd(x) = <Pcl(x + L) + Mn, Pci(xm) = nm (6) 

similarly to the standard density fluctuations in the Hal¬ 
dane scenario. Now, we can apply the standard repre- 













5 


sentation of the delta function 


<%(z)] 


E 

zeros of g 


\Vg(xj)\ 6{x 


Xj ) 


( 7 ) 


which, when combined with the previous ansatz, gives 
M 

p(x) = y f(x m )S(x - Xm) = 

m= 1 

M 

= ^Vci(x) y f(x m )S[ip c i(x) - -Km]. (8) 

m= 1 


The initial density formula is recovered by applying to 
the right-hand-side of the latter formula the delta func¬ 
tion transformation of Eq. ([7]). The next passage is then 
to Fourier transform the previous formula. This is rather 
different with respect to the Haldane scenario; in the lat¬ 
ter, one has to Fourier transform a standard Dirac comb. 
Here, the Dirac comb is in fact weighted by the clus¬ 
ter configuration f(x m ), which affects its Fourier compo¬ 
nents. Let us define the functional: 

M 

x[Vd{x)\ = y f(x m )S[ip c i(x) - Km] = (9) 

m= 1 
M 

= y S[ip c i(x) - 7tto] + y 5[ip c i{x) - 7rl], 
m=l eecl 


where the last sum is again performed on two-particle 
clusters. The first part can be Fourier transformed by 
considering the Poisson summation formula 

OO OO 

y S(x-nK)= y — e ~ i7rkx / K ( 10 ) 

n=— oo k=— oo 

thus leading to 



FIG. 3: Configurations allowed for the field ip(x). Panel (a): 
In the Luttinger liquid scenario, ip(x) can describe all pos¬ 
sible particle configurations and takes integer values at the 
position of each particle (red circles). The integrated particle 
density (thin orange line) follows the profile of ip. Panel (b): 
In the cluster Luttinger liquid scenario, p c i is constrained by 
the cluster structure, which assumes that a finite number of 
particles are packed into small clusters (blue circles). The in¬ 
tegrated particle density (thick orange line) does not follow 
the behavior of <p c z, but jumps whenever a two-particle cluster 
is encountered. 


where the new field ip c i (x)' represents fluctuations on top 
of the perfect cluster-crystalline solution of the problem, 
and er = M/N. In this way we obtain the final form of 
the mapping between the microscopic continuum density 
and the cluster fields 


N 


Pi - X) = m W V°i{x )<( y - f 


-iklfial(x) l _ 


(14) 


< k=—c 


’ —V tp c i (x)' < y a k e 

v k ——oo 


2 ik[tp c i (x)' — nncrx] 


The single particle operators can also be expressed by 
applying the same procedure of the Haldane formalism. 
In the bosonic case, we get: 


M oo 

lim Slipd(x) — Km] = _ e -*k<Pci(aO. m) 

M—>oo z ' 7T 

m=l k——oo 

The last important point is to re-absorb the last term 
in Eq. ([9]). By Fourier expanding all of its components, 
one gets renormalized c k coefficients in the previous ex¬ 
pression - not affecting its functional form. As such, we 
neglect these effects in the following, noticing that they 
are nevertheless expected to be very small when the num¬ 
ber of two-particle clusters is small, that is, M ~ N. We 
will then re-express the particle density as a function of 
the cluster operators as follows 

p(x) = ^Vipdix) { E (12) 

v k= —oo ) 


ip B (x) ~ yn — —Vip c i(x)'e ll3 ' 9oI ( x ' ) x 


y oi k e 


i2k[tp c i (x)' — irncrx] 


(15) 


< k =—c 


where the operator -& c i(x)' is the conjugate of ip' cl , satis¬ 
fying 


1 

7T 


Vpci(x)\ -ddly)' 


-iS(x - y ). 


(16) 


This can be verified by considering the role of the fac¬ 
tor /3; approximating the density operators with its non¬ 
oscillating part, we need the following commutation to 
hold 


where the numerical pre-factor has been introduced to 
compensate for the delta function renormalized coeffi¬ 
cients. We now proceed by rescaling the fields as follows: 

Pd{x) = —2p c i(x)' + 2nncrx, (13) 


-V<Pd(xy, e - i/3 ^ i( E = -6(x - y)e~ i ^^ v \ (17) 


which is satisfied if /3 = a 1 and Eq. (16) holds. The 
effective low-energy Hamiltonian obtained by introducing 
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the cluster fields in the microscopic Hamiltonian is thus 
a compactified boson theory 

H = T / dx + K ch (d x V cl ) 2 ] (18) 

which is gapless and conformal, with central charge c = 1 . 
All correlation functions of the microscopic operators can 
be evaluated using conventional techniques - see Ref . 14 
for a review. The theory is thus very similar to a TLL, 
with the exception that the mapping between micro¬ 
scopic and low-energy degrees of freedom presents re¬ 
markable differences. In the next subsection, we will 
illustrate the effects on correlations, while in Sec. m 
we will show numerically that the interpretation of the 
low-energy excitations differs considerably in a CLL with 
respect to the TLL. 


C. Correlation functions and structure factors in 
the cluster Luttinger liquid phase 


IV. PHASE DIAGRAM 

In this section, we investigate the phase transition from 
the TLL to the CLL phase by u sing n umerical simulations 
based on the DMRG algorithm 28 29 . As boundary effects 
are expected to be prominent in the cluster phase, we 
perform simulations using PBCs for systems up to L = 70 
sites, keeping up to 1400 states per block and with up to 
10 finite-size sweeps. The typical truncation error in the 
final sweep is of the order 2 • 10 ~ 6 for L < 60 and < 10 ~ 5 
for L = 70. Observables such as the local density rij 
deviate from the mean value n at most as rij—n < 3T0 -5 . 

In order to illustrate the generality of our findings, we 
investigate two relevant scenarios for the case r c = 2 , 
namely, different densities n = 3/7 and n = 4/10, which 
in the classical limit lead to different cluster densities as 
discussed in the previous sections (density regimes where 
dominant Umklapp terms appear have been investigated 
in Ref!®). In order to keep commensurability with the 
cluster structure, we considered sizes of (14,28,42, 56, 70) 
and (20,30,40,50,60,70) for the two cases, respectively. 


The mapping between the original operators and the 
emergent cluster fields in (141 and (15) allows us to make 


predictions on the scaling of both Green’s functions and 
density-density correlations. As an example, we focus 
here on the latter (the effects on the former have been 
discussed in Ref!®), which after bosonization take the 
form 


(p(z)p(0)) ~ n 2 


oi 02 cos(27m<ra;) 


x T ' 1 


+ ... , (19) 


with oi, a\ and 71 non-universal coefficients. This latter 
expression displays a radically different spatial modula¬ 
tion with respect to the standard Luttinger liquid sce¬ 
nario. In particular, assuming periodic boundary condi¬ 
tions (PBCs) peaks in the static structure factor 

S{q) = j e l(e ~ 3)q [( nerij) - n 2 ] (20) 

1,3 


A. Entanglement entropies 


In order to locate the transition point between the TLL 
and the CLL, as a first observable we monitor how the 
ground state entanglement properties change as a func¬ 
tion of V/t. We consider the bipartite von Neumann 
entanglement entropy 

S A = -Tr p A log p A , ( 22 ) 


where pa is the reduced density matrix of the sub-system 
A with respect to the rest of the chain. As both the TLL 
and the CLL are described at low-energies by a conformal 
field theor y, the entanglement entropy fulfills the follow¬ 
ing scalinjJ®® 


Sx(Z)=3 ln 


— sin(nl/L) 
7r 


+ C+ 0(1/1°), (23) 


are now not displayed at momenta associated with the 
particle density, as in TLLs. In particular, the lowest 
momentum peak is located at 


ki 


27t(1 — n) 

r c 


( 21 ) 


which for r c = 2 , gives k\ = 7 r(l — n). This is exactly 
the corresponding momentum peak of the classical clus¬ 
ter configuration, which is unrelated to the density of 
individual particles as in TLLs. We note that other sig¬ 
natures of departure from the TLL picture can also be 
found using level spectroscopy techniques 44 . 

In the next two sections, we provide numerical evidence 
for the existence of a CLL at intermediate couplings in 
Eq. Q, by monitoring the entanglement entropy, the 
spectral properties, and correlation functions. 


where L is the system size, l is the block length, C is a 
non-universal constant, and c is the central charge of the 
theory. Corrections of order 1/Z“, if any, are expected 
to exhibit no oscillations for the von Neumann entropy 
under PBCs. 

In order to extract the central charge of the system 
from finite-size simulations, we proceed in two-steps. 
First, we evaluate the finite-size value of the central 
charge, c(L), by fitting the entanglement entropy at a 
fixed size with the following function 


Sl(1) = < ^~K(l)+a 0 , k(1) = log 


sin( 7 r l/L) 


(24) 


where k(1) is the logarithm of the cord length. We keep 
l = 7 as the minimum block size considered in order 
to avoid possible corrections due to the breakdown of 
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FIG. 4: Entanglement entropy scaling as a function of the 
block length for n = 4/10. Panels (a)-(c): different interaction 
strengths V/t = 4.0,5.675,6.5. Data for different sizes are 
indicated by symbols, while lines are linear fits. Deep in the 
TLL phase, the finite-size central charge quickly approaches 
1. At the transition point, the value is 1.5. Deep in the 
CLL phase, finite-size-effects are strong - as evidenced by the 
bending of the entropy for different system sizes. Panel (d): 
L = 50, different interaction strengths. The black continuous 
and blue dot-dashed lines are guides for the eye at c = 1 and 
3/2, respectively. 



(b) 




FIG. 5: Entanglement entropy scaling as a function of the 
block length for n = 3/7. Panel (a): entropies at the tran¬ 
sition point V/t = 5.5. Data for different sizes are indicated 
by symbols, while the black line is a linear fit in good agree¬ 
ment with c = 3/2. Panel (b): L = 50, different interaction 
strengths. The black continuous and blue dot-dashed lines 
are guides for the eye at c = 1 and 3/2, respectively. 


Eq. (231 for blocks that are too small. Typical results are 
presented in Figs. [4] and [5] We find that a linear fit works 
quite well for all system sizes and in the entire parameter 
regimes we investigated. The finite-size central charge 
can be extracted with an accuracy of order 

A summary of all data is presented in Fig. |6ja), where 
we show the dependence of c(L) as a function of V/t for 
different system sizes and n = 4/10. The data strongly 
suggest that there is an intervening phase transition be¬ 
tween the two different c = 1 phases (TLL and CLL) at 
an intermediate value of V/t , evidenced by the bell-like 


FIG. 6: Panel (a): Finite-size central charge as estimated 
from a linear interpolation according to the Cardy-Calabrese 
formula. Blocks of size < 9 are neglected to avoid strong non- 
universal effects. Residuals are of order 10~ 4 in the worst 
cases, closer to the transition line and at large system sizes. 
The data evidence a phase transition around V/t ~ 5.63±0.02 
with a c = 3/2 central charge. Both liquid phases are com- 
pactified bosons, as expected. Panel (b): Finite-size scaling 
for the critical point V£ as a function of 1/L. Errors in the es¬ 
timates of V£, obtained via local interpolations of cl in panel 
(a), are of the size of the symbols. 


structure centered around V/t ~ 5.65. 

At the transition point, the central charge is compati¬ 
ble with c = 3/2, as shown in Fig. |4](b). This universal¬ 
ity class points toward the presence of a supersymmetric 
critical pointful, where the low-energy field theory is de¬ 
scribed by a combination of an compactified boson, and 
a real (Majorana) fermiorP^. Critical points of this kind 
have been recently discussed in different situations, such 
as constrained mod els com ing from explicitly supersym¬ 
metric HamiltonianJ- 31 32 -l and as effective boundary de¬ 
grees of freedom for topological phased. 

In order to extract the transition point more accu¬ 
rately, in Fig. |6](b) we plot as a function of the system 
size, the critical value of V£ at which the finite-size cen¬ 
tral charge reaches its maximum value, and then take 
the thermodynamic limit. We obtain an estimate of the 






























































critical point V c /t ~ 5.63 ± 0.02. 

B. Gaps and low-energy degrees of freedom 

In order to further deepen our understanding of the 
two phases and of the low-energy excitations in the vicin¬ 
ity of the transition point, we analyze the lower part of 
the spectrum of the Hamiltonian. In particular, we have 
targeted the single particle gap, defined as 

A sp (L) = 2 E n (L) - £jv-i(£) - E n+1 (L), (25) 

where En(L) is the ground state energy in a chain of 
length L with N particles. In addition, we have investi¬ 
gated the cluster gap A c ;, which, for the case n = 4/10, 
is defined as 

A cl (L) = 2 E n (L) - E n _ 2 (L) - E n+2 (L). (26) 

In the TLL phase, both gaps are expected to vanish, as 
they are linear combinations of different vertex operators. 
In the cluster phase, however, the picture is different, 
as discussed below. In the following, we first analyze 
the gaps in the classical limit t — 0 by considering 
the exact solution of Sec. [TT] and then present predic¬ 
tions of field theory and exact numerical results for t 0. 

a. Single-particle gap in a clustered state: Let us 
consider a system with size L = 10£, Na = 2£ and 
Ng = £, where £ labels the number of building blocks. 
The classical energy of the system is 

E n = VN b = Vi (27) 

and, upon doping, one gets 

E n+1 = V£ + 2V, E n _ 1 =Vi-V, (28) 

as the states with N + 1 and N — 1 particles cannot 
rearrange properly due to cluster constraints. This 
implies that the single particle gap is always open in the 
cluster phase, and in particular A sp (L) = V for every 
system size. 

b. Cluster gap in a clustered state: For the cluster 
gap the situation is remarkably different. Doping with a 
B cluster (two-particles) generates the following configu¬ 
ration in the classical limit 

E n+2 = V£ + 3V, E N — 2 = V£-3V (29) 

The latter is very much reminiscent of a classical crystal 
deformation - insertion and extraction of a single cluster 
take the same energy from the system. This leads to 
a vanishing cluster gap, as both contribution exactly 
cancel. 

From the field theory described in Subsec. |III B[ one 
can see that, while the cluster gap is nothing but a 


(a) (b) 



FIG. 7: Panel (a): Finite-size scaling of the cluster gap deep 
into the different phases. The gap scales to 0 in both CLL and 
TLL phases. Panel (b): data close to the transition point. 



FIG. 8: Finite-size scaling of the single particle gap deep into 
the different phases. The gap scales to 0 in the TLL phase, but 
takes a finite expectation value in the CLL phase, in agree¬ 
ment with the field theory predictions. 


combination of vertex operators in the cluster language 
and thus vanishes in the cluster phase, the single particle 
gap cannot be written in this way, thus implying that 
single particle excitations are never gapless within the 
low-energy description. This implies that there must be 
a phase transition between the CLL and TLL phases 
where the single particle gap A sp opens, while the 
cluster gap A c i remains largely unperturbed. 

The numerical results on the gap scaling fully confirm 
this picture. In Fig. [7ja), we show the cluster gap A c ; 
for some points representative of the TLL (V/t = 5) and 
cluster phase ( V/t = 6,7,8). In both phases, A d scales 
to 0 in the thermodynamic limit as a power law - as pre¬ 
dicted from the field theory analysis. Moreover, the scal¬ 
ing is also present close to the transition point [Fig. [7](b)] , 
where the gap vanishes approximately as oc 1/L. 

The behaviour of the single particle gap A sp in the two 
phases is illustrated in Fig. [8j Deep in the TLL phase, 
the gap scales to 0 (green diamonds and black pluses). 
However, once in the CLL phase, the gap is clearly finite, 
and extrapolates to larger values as V increases. We 
notice that in this latter region we have excluded from 
the fit to the data small system sizes with L = 10, which 
might display strong finite-size effects due to the presence 
of limited cluster structures at very small sizes. 
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FIG. 9: Finite-size scaling of the single particle gap in the 
vicinity of the transition point. Lines are best fits of the form 
ao + a.i(1/), with ao, a i and 02 constants; the caption 
indicates as a reference some of the extrapolated values. 



FIG. 10: Single particle gap extrapolated to the thermody¬ 
namic limit as a function of the interaction strength. The 
red line is a linear fit in the vicinity of the transition point, 
while the dashed line is a linear fit improved with a logarith¬ 
mic correction. Errors in the extrapolation are of order of the 
symbol-size (~ 0.01) except at V)t = 5.6,5.675, where the 
errors estimated by least-square methods are ~ 0.015. 


In Figs. [9| and [T0| we show how the single particle gap 
A sp scales in the vicinity of the transition point. In par¬ 
ticular, Fig. [TO] shows that the dependence of the gap 
on V/t is approximately linear (red dashed line), which 
is consistent with an emergent Ising field at the criti¬ 
cal point'®. The linear extrapolations locate the critical 
point around V/t ~ 5.605, in good agreement with the 
results for the scaling of the entanglement entropy (see 
Subsec. IVA). Moreover, we notice that the agreement 
improves considerably once additional logarithmic cor¬ 
rections are included in the fit. This could point towards 
the presence of additional corrections close to the critical 
point, which have already been observed in similar su¬ 
persymmetric scenarios 32 . An accurate understanding of 
the relation between the microscopic degrees of freedom 
and the emergent fields at the critical point might shed 
further light onto these finite-size corrections. 


C. Sound velocities 


At the critical point, the emergent supersymmetry im¬ 
plies that the scaling dimensions of the operators are 
fixed, and the sound velocities of both the boson and 
Ising (fermionic) mode are the same. Within our frame¬ 
work, we can identify the sound velocity of the bosonic 
mode vb using conventional conformal field theory tech¬ 
niques (see, e.g., Ref. ITT by monitoring the finite-size 
scaling of the cluster gap§" 


Ad/2 = 


27 tv b x b 
L 


(30) 


where Xb is the scaling dimension of the corresponding 
vertex operator. Following the scaling dimensions for the 
expected value of xb = 1/4A' at criticality, with K = 4 
being the Luttinger parameters!, one gets: 


4A C ,L 
V B = - 

7T 


(31) 


The estimate of the fermionic velocity is however not un¬ 
ambiguous, as our field theory does not allow to establish 
an exact mapping between the low-energy, continuum 
fields, and the lattice operators. As the gap in the Ising 
model approaches 1 in the strong coupling limit of the 
Ising model following the notations in Ref. 09, while our 
A sp approaches V and scales approximately as ~ (V—V c ) 
in the vicinity of the critical point, it seems plausible to 
assume that A sp = Ap in our context, where Ap is the 
gap in the Ising model corresponding to different parity 
sectors - thus related to the spin primary operator. This 
implies: 


A sp — 


2TTVf(S <7 + S(j ) 
L 


(32) 


where S a ,6a- = 1/16 are the anti-holomorphic and holo- 
morphic dimensions of the primary operator correspond¬ 
ing to the Ising spin at low energy. We thus get the 
following form for the fermionic sound velocity: 


v f = 


8A sp L 


(33) 


In Fig. 

Eq.pm 


we report our results for the velocities in 
in the vicinity of the critical point. The y- 


axis are rescaled by the value of the largest velocity in 
the L —00 limit to improve readability. The fits, indi¬ 
cated by the continuos lines, are 3-parameter fits of the 
form a,Q + ai* L° 2 , where the extrapolated value a 0 rep¬ 
resents the velocity in the thermodynamic limit. Typical 
errors in the extrapolated limit are of order 3/fp21. As 
can be seen in the upper panel, the velocities are equal 
(within numerical error) in the close vicinity of the tran¬ 
sition point, further supporting the supersymmetric na¬ 
ture of the critical point itself. We note that the next or¬ 
der corrections to the fermionic velocity are scaling with 
02 — —2, in full agreement with the conformal field the¬ 
ory predictiorf®. 
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1/L 


FIG. 11: Rescaled sound velocities Vp = v 13 /vm (where vm is 
the maximum of the two sound velocity in the L —> oo limit) 
of the fermionic and bosonic models are extracted from the 
low-energy spectrum following Eq. |31[ |32| 



FIG. 12 : The structure factor S(q) for V/t = 8, L = 56. For 
example, /x = — 3.2 1 (—5 1) for temperatures T = 1/t (8/t), 
such that IV ~ 24 is kept constant. 


In the following we focus on the static structure factor 
at finite temperature, in analogy to Eq. (201 


V. RESILIENCE OF CLUSTER FEATURES 
AGAINST TEMPERATURE 


S(q) = ^{(jii-n^rij -n)). (35) 

i,3 


In this section, we study the finite temperature prop¬ 
erties of the cluster Luttinger liquid state. Even though 
thermal fluctuations will suppress any long-range or 
quasi-long range order for the ID system at finite tem¬ 
perature, we will show that certain cluster features char¬ 
acterized by a peak in the static structure factor for cer¬ 
tain quasi-momenta can still survive for temperatures of 
the order of the interaction strength. To study the finite 
temperature properties of the system, we use a numeri¬ 
cally exact Quantum Monte Carlo (QMC) method with 
worm-type updated^ implemented in RefJ 51 : 

We choose the hard-core boson Hamiltonian we dis¬ 
cussed above with r c = 2: 


H = -t^2(b\b i+ 1 + h.c) + - 


Hn i 


. . !31! 

Here /x is the chemical potential, since our QMC simula¬ 
tions are performed in the grand canonical ensemble. We 
note that in Hamiltonian (34), t > 0 and the frustration 
induced by the next-nearest neighbour interactions only 
appears in diagonal terms. As a result, the system is 
free from the so-called ’’sign problem” irrespectively of 
the filling factor. In the QMC simulation, we choose the 
interaction strength V = 8t so that the ground state 
of the Hamiltonian corresponds to a cluster Luttinger 
liquid state. We focus on a chain with PBCs, and tune 
/i such that the average particle number is kept constant 
for the different temperatures studied. 


Here q is again the lattice quasi-momentum, n = N/L 
is the average density and (O) = Tr(Oe~^ H )/Z, with Z 
the partition function at temperature T = 1//3. For a 
chain with L = 56, N = 24 and PBCs the number of 
blocks of type A(B), N a(b), satisfy 

N a + 2 N b = N (36) 

3 N a + m B = L, (37) 


with N a = N b = 8. For these parameters, the cluster 
Luttinger liquid ground state exhibits a sharp peak in 
S(q) at a characteristic wave vector k\ = 7r(l — n) [from 
Eq. ©], which indicates cluster features. In Fig. [1^ 
we plot S(q) and tune /x such that the average number 
of particles is ( N ) ~ 24, for the different temperatures 
studied (T = 1/t, 8/t and 16/t). As expected, we find 
that the peaks, located at k\ = 47r/7, become broader 
with increasing T. However, they remain clearly visible 
up to temperatures of the order of V. This indicates that 
certain cluster features are comparatively robust against 
thermal fluctuations for sufficiently low temperatures. 


VI. ADIABATIC STATE PREPARATION OF A 
CLUSTER LUTTINGER LIQUID STATE 

If a system is prepared in the ground state of a certain 
initial Hamiltonian H 0 at tim e r = 0 and some parame¬ 
ter A(t) is adiabaticall} * 52 ! 53 ! tuned such that a different 
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Hamiltonian Hi is obtained at a final time r = T max , 
the adiabatic theorenP^ ensures that the original ground 
state will continuously evolve into the ground state of 
Hi. The total Hamiltonian H{r) can be written as 


H(t) = [1-X(t))H 0 + X(t)Hi, 


(38) 


where A (r) ranges from 0 to 1 as time runs from r = 0 

tO T = T max . 

As a prototypical example, let’s imagine to adi- 
abatically prepare the CLL state of a system with 
L = 14, N = 6 and r c = 2, that is the ground state of 

Hi = ~tiJ2i( b l b i+i + H - c -) + VJ2iJ2i=i n i n i+e for > 
e.g., V/ti = 6 (see Sec. IV). The initial Hamiltonian is 

H 0 = -~toJ2i( b l b i+i +H.c.) + wi th 

t 0 -C ti. Neglecti ng pe rturbatively small corrections 


oc ty/V (see Subsec. Ill A), the ground state of H 0 , under 


the assumption of PBC, is a mixed state of the 2-fold 
degenerate classical cluster configurations obtained from 
the possible permutations of two blocks of type A and 
two of type B. For L = 14, the allowed configurations 
are either A ABB or ABAB. 


Since the interaction part in H(t) is constant through¬ 
out the evolution, we can rewrite Eq. (38) as 


H(t) = H 0 + P(t)Hi, 


(39) 


where H 0 = V JV n%ni+u Hi = ~J2,( b l b i+i + 
H.c.) and /?(r) = t 0 (l - r/r max ) + ti(r/r max ). Notice 
that, as in Eq. 38 H( 0) = H 0 and H(r max ) = Hi. 


In Fig. 13 we show the static structure factor S(q), 
at r = r max , obtained by integrating the Schrodinger 
equation governed by the time-dependent Hamiltonian 
in Eq. (39) for different choices of V/ti, while keeping 


to/ti fixed. We compare results of the adiabatic pro¬ 
tocol with exact diagonalisation calculations of the final 
Hamiltonian Hi. 

Since both AABB and ABAB explicitly break the 
aforementioned ground state degeneracy, the corre¬ 
sponding peaks of S{q) have slightly different heights 
when compared to each other, despite being located 
at the same predicted critical momentum ki = 47t/ 7 
of Eq. (21). Independently of the choice of the initial 


state, we find that the best agreement between exact 
diagonalisation and state preparation is obtained in 


the strongly interacting limit [blue lines in Fig. 13 


where the CLL state is expected to be robust against 
quantum fluctuations (Sec. IV) and non-adiabatic effects. 


Experimentally, the protocol we propose here can be 
implemented by starting with particles trapped in a very 
deep lattice, and then by reducing the optical lattice 
depth with a down-ramp that ensures a linear growth in 
time of 0 (t). 


In order to take into account possible losses due to off- 
resonant light scattering induced by the optical lattice 



e—• V/ti = 2, ED 
-- AABB 

. ABAB 

«-♦ VI n = 4. ED 

-AABB 

. ABAB 

A—A V/ti = 6, ED 

- AABB 

. ABAB 


0 Tin lull 3jt/7 431/7 5n/7 6k 11 
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FIG. 13: Comparison of the adiabatic state preparation at the 
final time r ma x = 200/ii and starting from either ABAB (dot¬ 
ted lines) or AABB (dashed lines), with exact diagonalisation 
results (continuous lines plus symbols). Different colours dis¬ 
tinguish different Hamiltonian parameters V/t\. The ratio 
to/ti = 0.01 is fixed. Black squares indicate the classical 
limit prediction. 


and/or scattering of atoms with the background gas, we 
consider in addition a homogeneous and constant rate 7 
of particle losses. We assume that 7 takes into account 
also non-homogenous losses originated by excitations to 
Rydberg states, which might experience different Stark 
shifts with respect to ground state atoms. Three-body 
collisions can be safely neglected due to the hard-core 
assumption (tA ) 2 = 0. In the Born-Markov approxima¬ 
tion, the reduced system density matrix p(r) then evolves 
according to the master equation (H = 1 ) 


Hi 

p( T ) = -i(H ef£ (T)p(T) - p(t)hI s (t)) + 7 bjp(r)b 


3=1 


(4°) 

where H e g(r) = H(r) — (i/ 2)7 E J= i bjbj is an effective 
non-Hermitian Hamiltonian and H(r) is given by 
Eq. ([391). 


We approximate the exact evolution in Eq. (40) by 
Monte-Carlo trajector y sim ulations in the quantum jump 
appr oacl^-d 1 M) j n pjg 


14 


we show S{q) at r = r max , ob¬ 
tained by averaging over V tra j = 500 different system 
realisations. As expected, the best approximation of the 
non-dissipative case is given by the smaller choice of 7 . 
Nevertheless, the characteristic cluster peak in the static 
structure factor survives even for comparatively large loss 
rates, ensuring the stability of cluster features with re¬ 
spect to non-adiabatic effects in the protocol. 
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FIG. 14: Plot of S(q) for dissipative dynamics. The loss rates 
7 range over three order of magnitudes. Also the case of 
purely Hamiltonian dynamics (7 = 0) is included. The initial 
state is ABAB, r max = 200/ii and V/ti = 6 . Notice how 
cluster features survive at 7 = 10 ~ 3 ti. 

VII. CONCLUSION AND OUTLOOK 

In summary, we have provided an in-depth study of 
the emergence of cluster Luttinger liquid phases in ID 
models with soft-shoulder interactions. Starting from 
the exactly solvable classical limit, we have shown that 
different approaches such as strong coupling pertur¬ 
bation theory and constrained bosonization analysis 
predict the emergence of such a liquid phase of matter. 
In addition, we have analyzed numerically the full 
phase diagram for soft-shoulder potentials for different 
density regimes. Surprisingly, a critical point with 
central charge c = 3/2 is found to separate the CLL 
and TLL phases, which have both c = 1, indicating an 
emergent supersymmetry. While in the present work 
we have demonstrated that both the scaling of the 
entanglement entropy and of the gap provide evidence 
for such emergent behavior, and how the sound velocities 
of the emergent bosons and fermions are indeed equal 
within numerical accuracy, in future investigations 
it would be desirable to develop a held theoretical 
understanding of the emergent Ising held, and to provide 
a possible explanation of the stability of the critical 
point based on microscopic symmetries (as discussed in 
Ref. 1521) . We remark that c = 3/2 critical points have 
also been reported in bilinear-biquadratic spin-chains 
and multi-species Hubbard models in the presence of 
an explicit Z 2 symmetry related to binding-unbinding 


mechanisms.^^ 

In order to address questions of experimental ob¬ 
servability of the CLL phase in Rydberg-dressed gases, 
we have shown that finite temperature effects do not 
severely degrade the CLL signatures in correlation 
functions up to temperatures of the order of the 
interaction energy. We have further discussed how 
cluster states can adiabatically be prepared on realistic 
time-scales where the effect of decoherence is largely 
non-detrimental. Since experiments have already been 
performed with atomic Rydberg chains comprising up 
to 40 sites, these studies suggest that the realization of 
the desired system dynamics may be accessible within 
state-of-the-art technology. In particular, the departure 
from the conventional Luttinger liquid scenario could be 
benchmarked by monitoring correlation functions either 
via noise-correlation measurements or by extracting 
structure factors from single-site density measurements. 

An intriguing extension for the models discussed here 
is the two-dimensional scenario. There, the classical 
model is also exactly solvable, albeit the final Hamilto¬ 
nian looses the simple XY form derived in one-dimension. 
The possibility of emergent gauge fields and a violation of 
the Luttinger theorem in the fermionic case could lead to 
the stabilization of quantum spin liquid phases in a rela¬ 
tively simple model Hamiltonian, where cluster features 
play a prominent role. 
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The error is estimated by performing the fits with a dif- beyond the linear order which have non-universal nature 

ferent set of points, and keeping track of the variation in (see, e.g., the discussion in Ref. [49). 

the value of do- The term a\ * L a2 represents corrections 


